% Updates G for the Multisector GE model
%
% Jon Steinsson and Emi Nakamura, July 2006
%****************************************************

function GError = CalcGErrorMultiSectorCalvoPlus(Q,F,F2,G,alphaCalvo,...
    a,rp,rad,theta,Proba,ProbNgr,SectorWeights)

agridsize = size(a,1);
rpgridsize = size(rp,1);
radgridsize = size(rad,1);
Ssize = agridsize*rpgridsize*radgridsize;

NSectors = size(SectorWeights,1);

% Get weights as a function of S_t/P_t-1
%********************************************************

QaInd = repmat((1:agridsize)',rpgridsize*radgridsize,1);

QradInd = repmat((1:radgridsize)',[1 agridsize rpgridsize]);
QradInd = permute(QradInd,[2 3 1]);
QradInd = reshape(QradInd,[Ssize 1]);

OldInd = (1:Ssize)';

for ii = 1:NSectors
    [dummy,QrpInd] = ismember(int32(F(:,ii).F*10^6),int32(rp*10^6));
    clear dummy
    %QrpInd = match(int32(F(:,ii).F*10^6),int32(rp*10^6));
    FInd(ii).FInd = sub2ind([agridsize rpgridsize radgridsize],QaInd,QrpInd,QradInd);
    [FInd(ii).FInd iFInd] = sort(FInd(ii).FInd);
    OldFInd(ii).OldFInd = OldInd(iFInd);
    
    [dummy,QrpInd] = ismember(int32(F2(:,ii).F2*10^6),int32(rp*10^6));
    clear dummy
    %QrpInd = match(int32(F2(:,ii).F2*10^6),int32(rp*10^6));
    F2Ind(ii).F2Ind = sub2ind([agridsize rpgridsize radgridsize],QaInd,QrpInd,QradInd);
    [F2Ind(ii).F2Ind iF2Ind] = sort(F2Ind(ii).F2Ind);
    OldF2Ind(ii).OldF2Ind = OldInd(iF2Ind);
end
clear QaInd QrpInd1 QrpInd iFInd iF2ind OldInd

% Apply F
for ii = 1:NSectors
    Q_temp = accumarray(FInd(ii).FInd,Q(ii).Q(OldFInd(ii).OldFInd));
    Q_temp = [Q_temp; zeros(Ssize-size(Q_temp,1),1)];
    Q_temp2 = accumarray(F2Ind(ii).F2Ind,Q(ii).Q(OldF2Ind(ii).OldF2Ind));
    Q_temp2 = [Q_temp2; zeros(Ssize-size(Q_temp2,1),1)];
    Q(ii).Q = alphaCalvo(ii)*Q_temp + (1-alphaCalvo(ii))*Q_temp2;
end
clear Q_temp Q_temp2

% Apply Proba and ProbNgr
for ii = 1:NSectors
    Q(ii).Q = reshape(((reshape(Q(ii).Q,agridsize,rpgridsize*radgridsize)')*Proba(ii).Proba)',[Ssize 1]);
    Q(ii).Q = reshape(reshape(Q(ii).Q,[agridsize*rpgridsize radgridsize])*ProbNgr,[Ssize 1]);

end

% Get the prices
%**********************************************************

GInd = CalculateGInd(a,rp,rad,G);

for ii = 1:NSectors
    F(ii).F = F(ii).F(GInd);
    F2(ii).F2 = F2(ii).F2(GInd);
end

% Create Price Index
%***************************************************************

totWeight = zeros(radgridsize,1);
sumF = zeros(radgridsize,1);
for ii = 1:NSectors
    F(ii).F = exp(F(ii).F).^(1-theta);
    F2(ii).F2 = exp(F2(ii).F2).^(1-theta);
    F(ii).F = alphaCalvo(ii)*F(ii).F + (1-alphaCalvo(ii))*F2(ii).F2;
    Q(ii).Q = SectorWeights(ii)*Q(ii).Q;
    totWeight = totWeight ...
        + sum(reshape(Q(ii).Q,[agridsize*rpgridsize radgridsize]))';
    sumF = sumF + sum(reshape(F(ii).F.*Q(ii).Q,[agridsize*rpgridsize radgridsize]))';
end

% Deal with zero mass
W_zero = (totWeight < 10^(-10));
totWeight = totWeight + W_zero;
sumF_zero = (sumF < 10^(-10));
sumF = sumF + sumF_zero;

GError = log((sumF./totWeight).^(1/(1-theta)));
GError = (1-sumF_zero).*(1-W_zero).*GError ...
    + W_zero.*zeros(radgridsize,1) ...
    + (1-W_zero).*sumF_zero.*zeros(radgridsize,1);
